/* */

#include <math.h>
#include <vector>
#include "Interpolation2D.h"
#include "function.h"

void plot1()
{
    int m = 31, n = 61;
    float *xa = new float[m];
    float *ya = new float[n];
    float **za = new float*[m];
    for (int i = 0; i < m; i++)
        xa[i] = 1 + 0.1*i;
    for (int j = 0; j < n; j++)
        ya[j] = 2 + 0.1*j;        
    for (int i = 0; i < m; i++)
    {
        za[i] = new float[n];
        for (int j = 0; j < n; j++)
            za[i][j] = fsin(xa[i], ya[j]);
    }
    PolInt2D<double, float, 3> p = PolInt2D<double, float, 3>();
    p.construct(xa, ya, za, m, n);
    for (float y = 0.0; y < 10.05; y += 0.05)
    {
        for (float x = 0.0; x < 5.05; x += 0.05)
            printf("%f,", p.evaluate(x, y));
        printf("\n");
    }
    for (int i = 0; i < m; i++)
        delete[] za[i];
    delete[] xa, ya, za;
}

void plot2()
{
    int m = 7, n = 7;
    float xa[] = {0.00117378,0.00472397,0.0171859,0.0409838,0.0780596,0.130081,0.198535};
    float ya[] = {0.365616,1.893040,3.786090,7.572170,15.144300,30.288700,60.577400};
    float za[7][7] = {
        0.0217822,0.0253224,0.0288237,0.0346827,0.0448323,0.0636086,0.100366,
        0.0233179,0.0268545,0.0303556,0.0362159,0.0463659,0.0651426,0.101902,
        0.0296429,0.0331470,0.0366371,0.0425000,0.0526603,0.0714467,0.108208,
        0.0402311,0.0440292,0.0477457,0.0538394,0.0641187,0.0829203,0.119654,
        0.0511250,0.0554077,0.0595859,0.0662932,0.0771901,0.0963434,0.133061,
        0.0625876,0.0673198,0.0719785,0.0794046,0.0910973,0.110757,0.147656,
        0.0748282,0.0800098,0.0851434,0.0933663,0.106111,0.126669,0.163872};
    float *zp[7];
    for (int i = 0; i < 7; i++)
        zp[i] = &za[i][0];
    SplInt2D<double, float> p = SplInt2D<double, float>();
    p.construct(xa, ya, zp, m, n);
    // for (float y = 0.0; y < 70.1; y += 0.1)
    // {
    //     for (float x = 0.0; x < 0.21; x += 0.01)
    //         printf("%f,", p.evaluate(x, y));
    //     printf("\n");
    // }
    std::vector<float> xv = 
        {0.0001,0.0003,0.0008,0.0022,0.0060,0.0167,0.0464,0.1292,0.3594,1.0000};
    std::vector<float> yv = 
        {0.1000,0.2154,0.4642,1.0000,2.1544,4.6416,10.0000,21.5443,46.4159,100.0000};
    for (float y : yv)
    {
        for (float x : xv)
            printf("%f,", p.evaluate(x, y));
        printf("\n");
    }

}

int main()
{
    plot2();
}